rm(list = ls())
require(stargazer)
require('clusterSEs')
require(mfx)


#Load in some Data
data2=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type2clean.csv", header=TRUE)
data3=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type3clean.csv", header=TRUE)
data4=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type4clean.csv", header=TRUE)
data5=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type5clean.csv", header=TRUE)

data16=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type16clean.csv", header=TRUE)
data17=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type17clean.csv", header=TRUE)
data18=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type18clean.csv", header=TRUE)
data19=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type19clean.csv", header=TRUE)

dataframe16=data.frame(data16)
dataframe16=subset(dataframe16,uid>1600)
dataframe17=data.frame(data17)
dataframe17=subset(dataframe17,uid>1700)
dataframe18=data.frame(data18)
dataframe18=subset(dataframe18,uid>1800)
dataframe19=data.frame(data19)
dataframe19=subset(dataframe19,uid>1900)

dataframe2=rbind(data.frame(data2),dataframe16)
dataframe3=rbind(data.frame(data3),dataframe17)
dataframe4=rbind(data.frame(data4),dataframe18)
dataframe5=rbind(data.frame(data5),dataframe19)

#Properly assign parameters to different experiment types
dataframe2$version=1
dataframe2$Bcorrinc=55
dataframe2$Bfailinc=40


dataframe3$version=2
dataframe3$Bcorrinc=52
dataframe3$Bfailinc=40


dataframe4$version=3
dataframe4$Bcorrinc=52
dataframe4$Bfailinc=30

dataframe5$version=4
dataframe5$Bcorrinc=55
dataframe5$Bfailinc=30

#Combine data and assign useful dummies
dataframe=rbind(dataframe2,dataframe3,dataframe4,dataframe5)
dataframe$Bresponse=(dataframe$response==19)
dataframe$Bstate=(dataframe$state==4)
dataframe$Cpresent=(dataframe$qtype==4)

#Show N for each experiment type
#**
length(unique(dataframe2$uid))
length(unique(dataframe3$uid))
length(unique(dataframe4$uid))
length(unique(dataframe5$uid))
#**

#start making table

Bperc49C=vector('numeric')
Bperc49noC=vector('numeric')
diffp49=vector('numeric')

Bperc51C=vector('numeric')
Bperc51noC=vector('numeric')
diffp51=vector('numeric')
pecentsig=vector('numeric')
subjectperc=vector('numeric')

for(i in 1:4){
#Aggregate tests
versionframe=subset(dataframe,version==i)
dataframe49=subset(versionframe,state==3)
Bperc49C[i]=sum(dataframe49$Cpresent*dataframe49$Bresponse)/sum(dataframe49$Cpresent)
Bperc49noC[i]=sum((1-dataframe49$Cpresent)*dataframe49$Bresponse)/sum(1-dataframe49$Cpresent)
model49=logitmfx(Bresponse~Cpresent,data=dataframe49, robust=TRUE, clustervar1 = "uid") 
stargazer(model49$mfxest)
stargazer(model49$fit)
diffp49[i]=model49$mfxest[4]

dataframe51=subset(versionframe,state==4)
Bperc51C[i]=sum(dataframe51$Cpresent*dataframe51$Bresponse)/sum(dataframe51$Cpresent)
Bperc51noC[i]=sum((1-dataframe51$Cpresent)*dataframe51$Bresponse)/sum(1-dataframe51$Cpresent)
model51=logitmfx(Bresponse~Cpresent,data=dataframe51, robust=TRUE, clustervar1 = "uid") 
stargazer(model51$mfxest)
stargazer(model51$fit)
diffp51[i]=model51$mfxest[4]

treatmentuidlist=unique(dataframe51$uid)
indiprob=vector('numeric')
for(j in 1:length(treatmentuidlist)){
indiframe=subset(dataframe51,uid==treatmentuidlist[j])
model51indi=logitmfx(Bresponse~Cpresent,data=indiframe, robust=TRUE) 
#note additive term insures wrong directions gets invalid pvalue
indiprob[j]=model51indi$mfx[4]+(model51indi$mfx[1]<=0)
}
subjectperc[i]=100*sum(indiprob<0.05)/length(indiprob)
print(i)
}

i=5
versionframe=dataframe

dataframe49=subset(versionframe,state==3)
Bperc49C[i]=sum(dataframe49$Cpresent*dataframe49$Bresponse)/sum(dataframe49$Cpresent)
Bperc49noC[i]=sum((1-dataframe49$Cpresent)*dataframe49$Bresponse)/sum(1-dataframe49$Cpresent)
model49=logitmfx(Bresponse~Cpresent,data=dataframe49, robust=TRUE, clustervar1 = "uid") 
stargazer(model49$mfxest)
stargazer(model49$fit)
diffp49[i]=model49$mfxest[4]

dataframe51=subset(versionframe,state==4)
Bperc51C[i]=sum(dataframe51$Cpresent*dataframe51$Bresponse)/sum(dataframe51$Cpresent)
Bperc51noC[i]=sum((1-dataframe51$Cpresent)*dataframe51$Bresponse)/sum(1-dataframe51$Cpresent)
model51=logitmfx(Bresponse~Cpresent,data=dataframe51, robust=TRUE, clustervar1 = "uid") 
stargazer(model51$mfxest)
stargazer(model51$fit)
diffp51[i]=model51$mfxest[4]

#Individual Tests
treatmentuidlist=unique(dataframe51$uid)
indiprob=vector('numeric')
for(j in 1:length(treatmentuidlist)){
indiframe=subset(dataframe51,uid==treatmentuidlist[j])
model51indi=logitmfx(Bresponse~Cpresent,data=indiframe, robust=TRUE) 
#note additive term ensures wrong directions gets invalid pvalue
indiprob[j]=model51indi$mfx[4]+(model51indi$mfx[1]<=0)
}
subjectperc[i]=100*sum(indiprob<0.05)/length(indiprob)


rowlabels=data.frame(Treatment=c(1,2,3,4,"total"))
tablestate49=data.frame(PofBnoC=Bperc49noC,PofBC=Bperc49C,Prob=diffp49)
tablestate51=data.frame(PofBnoC=Bperc51noC,PofBC=Bperc51C,Prob=diffp51)


fulltable=cbind(rowlabels,tablestate49,tablestate51)
fulltable$subject_fraction=subjectperc

fulltable

